Skip to content

InnerLayer.SLAYER - FEATURE! - Build P_perp from the perpendicular particle diffusivity by default and add an ion heat diffusivity profile - #501

Open
ebursch wants to merge 15 commits into
developfrom
feature/slayer-pperp-models
Open

ebursch wants to merge 15 commits into
developfrom
feature/slayer-pperp-models

Conversation

@ebursch

@ebursch ebursch commented Oct 5, 2026 •

Copy link
Copy Markdown
Collaborator

Release note

  • Audience: users
  • Numerical impact: SLAYER layer parameters and growth rates move on every deck. P_⊥ is now built by default from the perpendicular particle diffusivity D_⊥, which on diiid_slayer_n1 is 0.22–0.33 of the old electron-χ value, so γ moves by +31 % (2/1), +18 % (3/1) and +24 % (4/1) while ω moves by ≤ 0.4 %. The new analytic χ profiles of the DIII-D-like kinetic file alone move γ by ≤ 1.5 %. The layer-overlap cap of diiid_truncation_n1 is unchanged (its widths move 2 %). (harness @ 4ce6bfe)
  • Migration: in [SLAYER], rename chi_perp to chi_perp_e (also the chi_perp keyword of build_slayer_inputs and slayer_parameters). The new default P_perp_model = "D_perp" needs the ion diffusivity χ⊥,i: add a chi_i profile to the kinetic file or set chi_perp_i, or set P_perp_model = "chi_perp_e" to keep the previous P_⊥. tau_E is the whole-plasma confinement time.

SLAYER's perpendicular Prandtl number P_⊥, which enters the inner-layer Δ(Q) with P_φ, can now be built several ways through [SLAYER] P_perp_model, and by default comes from the perpendicular particle diffusivity of Fitzpatrick (2022), Eq. 16. The critical Δ keeps the electron χ⊥,e in every mode, so the two are no longer forced to be the same diffusivity. The kinetic HDF5 schema gains an optional chi_i profile.

P_perp_model P_⊥ Reference
D_perp (default) τ_R D_⊥/r_s², with D_⊥ = c_β² η/μ₀ + (2/3)(1 − c_β²)(n T_e′ χ⊥,e + n T_i′ χ⊥,i)/p′ Fitzpatrick, Phys. Plasmas 29, 032507 (2022), Eqs. 4–6, 16
chi_perp_e τ_R χ⊥,e/r_s² (the only option before; Fitzpatrick's TJ convention) —
chi_perp_i τ_R χ⊥,i/r_s² —
P_phi P_⊥ = P_φ Fitzpatrick, Phys. Plasmas 30, 092512 (2023), Eq. 143
c_beta P_⊥ = c_β², the χ⊥ → 0 limit of D_perp J.-K. Park, Phys. Plasmas 29, 072506 (2022), Eq. 9 with K = 0
tau_E P_⊥ = P_φ = τ_R/τ_E(r_s), with τ_E(r) = r²/χ_E and a uniform χ_E = (3/16)a²/τ̄_E Fitzpatrick (2023), Eqs. 37–38, 81–83, 131

What changed

  • P_⊥ models in slayer_parameters, selected by [SLAYER] P_perp_model and passed through TearingProblem and the layer-overlap scan, so the scan and the tearing stage see one plasma.
  • D_perp default. A run that needs χ⊥,i and has neither a chi_i profile nor chi_perp_i stops with a message naming both. slayer_parameters keeps chi_perp_e as its own default, since it has no profile gradients.
  • tau_E takes the whole-plasma τ̄_E, converted to a uniform χ_E with κ_e = 3/16 (the T_e ∝ 1 − r̂² value of Eq. 38) and the boundary minor radius from the flux label √(2ψ_t/B₀).
  • Kinetic file. TkMkr_D3Dlike_Hmode_kinetic.h5 carries analytic χ profiles (below) in place of profiles derived from a measured discharge, and gains chi_i; n, T and ω_E are byte-identical to develop's. The forms and their sources are in the file's provenance attribute and in the stability docs.
  • Truncation deck pins P_perp_model = "chi_perp_e": under D_perp the DIII-D-like layers are about 30 % narrower and do not overlap inside ψ = 1, so the pin keeps the layer-overlap cap exercised.
  • Harness case diiid_slayer_n1_chi_perp_e keeps a baseline for the previous model.
  • Docs: the P_⊥ models and the analytic χ in stability.md, chi_i in the kinetic schema, the Fitzpatrick 2022 and 2023 papers in docs/resources, and the Park citation in the SLAYER file headers corrected to Phys. Plasmas 29, 072506 (2022).
  • An old deck's chi_perp is rejected with a hint that it is now chi_perp_e.
chi_e   = 2.75 + ψ + 3 exp(−25 (ψ − 0.85)²)        [m²/s, ψ = ψ_N]
chi_i   = 1/2 + (3/2) ψ²⁰ + (5/4) exp(−25 (ψ − 0.7)²)
chi_phi = 3/2 − (5/4) (1 − ψ)⁴

These are not fits to a discharge. Their magnitudes are of the same order as the TRANSP power-balance χ_e and χ_i of DIII-D H-modes in Hong et al., Phys. Plasmas 30, 072512 (2023), Fig. 5; χ_e keeps the factor-2 raise that makes the default deck exercise the overlap cut; χ_φ tracks χ_i (Prandtl number of order one, rising outward) as reported for DIII-D by Yoshida et al., Nucl. Fusion 52, 123005 (2012). They do not resolve the pedestal transport barrier, where power-balance χ falls to about 0.1–0.3 m²/s (Callen et al., Nucl. Fusion 50, 064004 (2010)).

Regression report

regress --refs 9fb32c9a6,4ce6bfea5 --force on feynman, all 20 cases (develop against this branch with develop merged). Tests and the docs build pass at the head.

case result
diiid_slayer_n1 P_⊥ ×0.22–0.33 at all six surfaces; γ 222.6 → 291.0 (2/1), −288.5 → −339.5 (3/1), −1124 → −1395 s⁻¹ (4/1); ω ≤ 0.4 %; no root flag changes
diiid_slayer_n1_chi_perp_e (new) no develop run (develop rejects the new key); at this head its γ is within 0.8 %, 0.9 % and 1.5 % of develop's diiid_slayer_n1, the effect of the analytic χ alone
diiid_truncation_n1 cap, psilim 0.99498, qlim 6.2 and surfaces unchanged; Eq. 100 widths 2.05 %
every other case unchanged
Full regression report
Regression Report: diiid_slayer_n1
===========================================================================================
Ref 1: 9fb32c9a6  @ 9fb32c9a6 (2026-10-20)
       env: julia 1.11.9, x86_64-linux-gnu, manifest 8e24ca66 (pinned), 8 threads/8 BLAS, jit build
Ref 2: 4ce6bfea5  @ 4ce6bfea5 (2026-10-20)
       env: julia 1.11.9, x86_64-linux-gnu, manifest 8e24ca66 (pinned), 8 threads/8 BLAS, jit build
-------------------------------------------------------------------------------------------
Quantity                            9fb32c9a6  4ce6bfea5  Diff                Status
-------------------------------------------------------------------------------------------
SLAYER surface indices              [6 elem]   [6 elem]   0.0e+00             OK
SLAYER poloidal m                   [6 elem]   [6 elem]   0.0e+00             OK
SLAYER toroidal n                   [6 elem]   [6 elem]   0.0e+00             OK
SLAYER minor radius rs              [6 elem]   [6 elem]   0.0e+00             OK
SLAYER r-based shear                [6 elem]   [6 elem]   0.0e+00             OK
SLAYER Lundquist S                  [6 elem]   [6 elem]   0.0e+00             OK
SLAYER D_norm                       [6 elem]   [6 elem]   0.0e+00             OK
SLAYER P_perp                       [6 elem]   [6 elem]   5.589e+01 (77.66%)  ** CHANGED **
SLAYER tauk                         [6 elem]   [6 elem]   0.0e+00             OK
SLAYER iota_e                       [6 elem]   [6 elem]   0.0e+00             OK
SLAYER Q_e                          [6 elem]   [6 elem]   0.0e+00             OK
SLAYER Q_i                          [6 elem]   [6 elem]   0.0e+00             OK
SLAYER Q_root [2/1,3/1,4/1]         [3 elem]   [3 elem]   6.787e-03 (1.09%)   ** CHANGED **
SLAYER ω_Hz [2/1,3/1,4/1]           [3 elem]   [3 elem]   4.435e+01 (0.40%)   ** CHANGED **
SLAYER γ_Hz [2/1,3/1,4/1]           [3 elem]   [3 elem]   2.713e+02 (24.14%)  ** CHANGED **
SLAYER no_root flags [2/1,3/1,4/1]  [3 elem]   [3 elem]   0.0e+00             OK
SLAYER enabled flag                 1          1          0.0e+00             OK
Runtime (s)                         298.2s     262.5s                         --
===========================================================================================
Summary: 4 changed, 13 unchanged

Regression Report: diiid_truncation_n1
=====================================================================================================
Ref 1: 9fb32c9a6  @ 9fb32c9a6 (2026-10-20)
       env: julia 1.11.9, x86_64-linux-gnu, manifest 8e24ca66 (pinned), 8 threads/8 BLAS, jit build
Ref 2: 4ce6bfea5  @ 4ce6bfea5 (2026-10-20)
       env: julia 1.11.9, x86_64-linux-gnu, manifest 8e24ca66 (pinned), 8 threads/8 BLAS, jit build
-----------------------------------------------------------------------------------------------------
Quantity                               9fb32c9a6      4ce6bfea5      Diff               Status
-----------------------------------------------------------------------------------------------------
layer-overlap cap applied              1              1              0.0e+00            OK
layer-overlap psilim                   9.985081e-01   9.985081e-01   0.0e+00            OK
first overlapping surface index        7              7              0.0e+00            OK
diffusive-resistive layer widths (Δψ)  [7 elem]       [7 elem]       1.320e-04 (2.05%)  ** CHANGED **
qlim                                   6.200000e+00   6.200000e+00   0.0e+00            OK
psilim                                 9.949767e-01   9.949767e-01   0.0e+00            OK
mpert                                  37             37             0.0e+00            OK
# singular surfaces                    5              5              0.0e+00            OK
singular q values                      [5 elem]       [5 elem]       0.0e+00            OK
singular psi locations                 [5 elem]       [5 elem]       0.0e+00            OK
plasma energy Re(ep[1])                -1.372382e+00  -1.372382e+00  0.0e+00            OK
total energy Re(et[1])                 8.005423e-01   8.005423e-01   0.0e+00            OK
Runtime (s)                            227.1s         199.9s                            --
=====================================================================================================
Summary: 1 changed, 11 unchanged

diiid_slayer_n1_chi_perp_e fails on develop, which does not know P_perp_model, so the harness exits with one failed run; gal_resistive_pe extracts nothing on either ref, as on develop. Every other case is unchanged.

Notes for reviewers

  • Taken over from @ebursch while he is on leave. His commits are kept; the develop merges, the review changes, the D_perp default, the analytic χ and the docs are on top.
  • Why D_⊥ for both the layer and Eq. 100. Fitzpatrick (2025) derives the diffusive-resistive width (Eq. 100) from the layer equations of Fitzpatrick (2022), where P_⊥ = τ_R/τ_⊥ with τ_⊥ = r²/D_⊥ multiplies the diffusion of the density variable (2022, Eq. 29). The 2025 paper writes χ_E in that slot and sets it to 1 m²/s, relying on the 2022 finding that D_⊥ is comparable to the energy diffusivities. Eq. 16 makes that precise; on these profiles it gives D_⊥ = 0.22–0.32 χ_e at the outer rational surfaces, because the weights n T_s′/p′ sum to less than one where density gradients carry part of p′.
  • Every P_⊥ model is μ₀χ/η, so r_s cancels except through η's trapped fraction. Only tau_E carries a radial label, through a; it uses the flux label whatever rs_method is, so the pending radial-label default change (InnerLayer.SLAYER - API - Choose one default radial label for the Fitzpatrick layer formalism #417) does not touch these models.
  • κ_e is fixed at 3/16 rather than computed from the profile, because Eq. 38 needs dT_e/dr at the edge, which a pedestal makes ill-conditioned (0.045–0.12 on this profile depending on how the edge slope is taken). The docstring states the order-2–3 shape uncertainty, and that in H-mode the χ_E implied by τ̄_E falls below the core power-balance χ.
  • n and T in the example file are unchanged: they are the long-standing DIII-D-like H-mode profiles, rough fits that every pinned harness quantity depends on.
  • Overlap with InnerLayer.SLAYER - DOCS - 🌱 Describe chi_tor as the toroidal momentum diffusivity, not a heat diffusivity #471, which edits the chi_tor docstrings in the same files: whichever merges second needs a small conflict resolution.

NO MERGE WITHOUT HUMAN REVIEW. This PR needs a third-party human reviewer's approval of the head that lands.

🤖 Generated with Claude Code

ebursch and others added 5 commits October 5, 2026 10:28
… kinetic HDF5

Add an optional chi_i (ion perpendicular heat diffusivity) dataset to the
kinetic-profile reader and writer, alongside chi_e and chi_phi.

Ship chi_i in the DIII-D-like example kinetic file, built like its chi_e:
AOT CHI_I of DIII-D 147131 at 2300 ms, mapped rho_N -> psi_N with EFIT01
RHOVN, the spurious zero at psi_N = 1 dropped, Gaussian-smoothed with
sigma = 12 of 201 rho_N points (the width that best reproduces the
example chi_e from AOT CHI_E). Range 0.46-2.1 m^2/s.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…i_perp to chi_perp_e

New [SLAYER] P_perp_model selects how the perpendicular Prandtl number is built:
- chi_perp_e (default, unchanged numbers): tau_R*chi_perp_e/r_s^2
- chi_perp_i: tau_R*chi_perp_i/r_s^2 (kinetic-file chi_i, scalar fallback)
- P_phi: P_perp = P_tor (Fitzpatrick, Phys. Plasmas 30, 092512 (2023), Eq. 143)
- D_perp: tau_perp = r_s^2/D_perp (Fitzpatrick, Phys. Plasmas 29, 032507 (2022))
- c_beta: C^2 = P_perp with C = c_beta (Park et al., Phys. Plasmas 29, 122505 (2022))

Breaking: the [SLAYER] chi_perp key is renamed chi_perp_e. chi_perp_e still sets
the critical-Delta offset in every mode. Regression harness diiid_slayer_n1
develop vs local: 17 unchanged.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…= tau_R/tau_E)

A single [SLAYER] tau_E [s] for the whole plasma sets P_perp = P_tor = tau_R/tau_E
at every surface (Fitzpatrick, Phys. Plasmas 30, 092512 (2023), Eq. 131).
validate() requires a positive tau_E when this model is selected.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…. 16

P_perp_model = D_perp now builds D_perp per surface from Eq. 16 of Fitzpatrick,
Phys. Plasmas 29, 032507 (2022), with eta_perp = eta_par:
  D_perp = c_beta^2 eta/mu0 + (2/3)(1 - c_beta^2)[eta_e tau_e/(1+eta_e) chi_perp_e
                                                 + eta_i tau_i/(1+eta_i) chi_perp_i]
using the n, T_e, T_i log-gradients at the surface (Eqs. 4-6, written flat-density
safe). The scalar D_perp input is removed.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…e DIII-D-like SLAYER example

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
@github-actions github-actions Bot added changed-results Results move or an interface breaks - read before upgrading feature New capability labels Oct 5, 2026
@ebursch ebursch self-assigned this Oct 5, 2026
@ebursch

ebursch commented Oct 5, 2026

Copy link
Copy Markdown
Collaborator Author

Effect of each P_perp_model on the DIII-D-like SLAYER example

Setup: examples/DIIID-like_SLAYER_example at 27d0744, run once per model with only P_perp_model changed. The run is uncoupled, uses AMR and is validity-gated. χ⊥,e, χ⊥,i and χ_φ come from the chi_e, chi_i and chi_phi profiles in TkMkr_D3Dlike_Hmode_kinetic.h5. For tau_E I used τ_E = 0.1 s, a typical DIII-D H-mode value, which is not derived from this equilibrium. All three surfaces found a root in every run.

P_perp_model P_⊥ (2/1, 3/1, 4/1) P_φ (2/1, 3/1, 4/1) γ [Hz] (2/1, 3/1, 4/1) Δγ vs default (2/1, 3/1, 4/1) ω [Hz] 2/1
chi_perp_e 50.3, 44, 35.7 34.7, 24.7, 19 192, -251, -1012 — -10527
chi_perp_i 29.9, 26.4, 12.5 34.7, 24.7, 19 213, -274, -1217 +11 %, -9 %, -20 % -10523
P_phi 34.7, 24.7, 19 34.7, 24.7, 19 206, -277, -1123 +8 %, -10 %, -11 % -10525
D_perp 16.4, 17.5, 10.9 34.7, 24.7, 19 242, -296, -1249 +26 %, -18 %, -23 % -10516
c_beta 0.0133, 0.0088, 0.00645 34.7, 24.7, 19 313, -588, -1669 +63 %, -134 %, -65 % -10133
tau_E 34.5, 30.8, 27.6 34.5, 30.8, 27.6 206, -264, -1034 +8 %, -5 %, -2 % -10525

The Δγ column is (γ − γ_default)/|γ_default|, so + means more unstable. For scale, the run-to-run AMR noise on γ is about 0.01 %.

Takeaways

  • Only c_beta changes the answer qualitatively. Its P_⊥ = c_β² ≈ 0.01, three to four orders of magnitude below every other model. The 2/1 drive rises about 60 %, the stable 3/1 and 4/1 roots damp 65–135 % harder, and ω on the 2/1 shifts about 4 %. Every other model shifts ω by less than 0.1 %.
  • The χ-based models move γ by roughly 10–25 %.
    • chi_perp_i uses the example's χ⊥,i, which is about 0.4–0.6 of χ⊥,e at these surfaces. That halves P_⊥ and moves γ by about 10 % at the 2/1 and 3/1 and 20 % at the 4/1.
    • D_perp gives the smallest P_⊥ of the transport-based models, so it moves γ the furthest, by 18–26 %.
    • Across these models, a lower P_⊥ makes the 2/1 more unstable and the 3/1 and 4/1 more stable.
  • P_phi and tau_E land close together. With τ_E = 0.1 s, τ_R/τ_E happens to be close to the χ_φ-based P_φ on the 2/1. tau_E also overrides P_φ, so its 3/1 and 4/1 sit nearer the default. γ scales with the chosen τ_E, so these numbers only illustrate the model.
  • chi_perp_e (the default) reproduces develop exactly; see the regression report in the PR description.

The unit tests (Pkg.test() at 27d0744 on feynman) pass.


🤖 Generated with Claude Code

@ebursch

ebursch commented Oct 5, 2026 •

Copy link
Copy Markdown
Collaborator Author

Example kinetic-file χ profiles vs the DIII-D 147131 AOT source

chi profiles: example .h5 vs AOT

Each panel compares TkMkr_D3Dlike_Hmode_kinetic.h5 (red dashed) with the AOT profile it was built from (blue): DIII-D 147131 at 2300 ms, mapped ρ_N → ψ_N with EFIT01 RHOVN. The grey lines mark the 2/1, 3/1 and 4/1 surfaces used in the SLAYER example. The bottom-right panel overlays the three .h5 profiles on one axis.

Dataset Status Source
chi_e (χ⊥,e) existing, unchanged AOT CHI_E, smoothed
chi_phi (χ_φ) existing, unchanged AOT ANG_MOM_DIFF, smoothed
chi_i (χ⊥,i) added in this PR AOT CHI_I; spurious zero at ψ_N = 1 dropped; Gaussian-smoothed with σ = 12 of 201 ρ_N points

At the SLAYER surfaces (new chi_i):

Surface ψ_N chi_i (.h5) AOT CHI_I .h5 vs AOT chi_i / chi_e
2/1 0.518 1.13 1.06 +7 % 0.59
3/1 0.770 1.80 2.10 −15 % 0.60
4/1 0.893 1.15 0.82 +39 % 0.35

Notes

  • Smoothing width: σ = 12 is the width that best rebuilds the existing chi_e from AOT CHI_E, about 10 % RMS in log over the profile.
  • Edge treatment differs: the existing chi_e and chi_phi also smooth over the pedestal well, the edge spike in CHI_E and the negative ANG_MOM_DIFF beyond ψ_N ≈ 0.94, so they stay positive and smooth to the edge. The file doesn't record how those edge treatments were done. The new chi_i uses only the Gaussian filter, so it follows the AOT well near ψ_N ≈ 0.93 and the rise to about 2 m²/s at the edge more closely.
  • 4/1 surface: the 4/1 sits on the inner flank of that well, where the smoothing lifts chi_i 39 % above AOT. A narrower σ would track AOT more closely there, at the cost of keeping more of the AOT structure elsewhere.

🤖 Generated with Claude Code

ebursch and others added 3 commits October 5, 2026 13:14
… and describe c_beta as the no-transport limit of D_perp

The SLAYER file headers and the c_beta P_perp_model cited "Park et al.,
Phys. Plasmas 29, 122505 (2022)", which does not exist. The paper is
J.-K. Park, Phys. Plasmas 29, 072506 (2022), doi:10.1063/5.0093079, as
docs/src/citations.md already lists. Its Eq. 10 defines
C² = c_β² + (1 − c_β²)K; C = c_β drops the thermal-conduction K, which
makes P_perp = c_β² the χ⊥ → 0 limit of the D_perp model.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
@ebursch
ebursch marked this pull request as ready for review October 5, 2026 20:11
@ebursch
ebursch requested a review from d-burg October 5, 2026 20:11
@ebursch

ebursch commented Oct 5, 2026 •

Copy link
Copy Markdown
Collaborator Author

@logan-nc @d-burg This adds a variety of options for perpendicular diffusion inputs. It doesn't change the default yet, so selecting the chi_perp_e is the same as what currently exists. You can see in the comparison that the only really different one is the c_beta which is good news (and encouragement to have this increased fidelity vs the Fortran slayer implementation). I think that we should move to full D_perp as default since there is not a literature justification for chi_perp_e ~ D_perp and if people already need chi_phi and chi_perp_e inputs than its not too crazy to also have them input chi_perp_i. It's encouraging that the tau_E choice is right order of magnitude (using tau_E = chi_phi and P_perp = P_phi) so that is also a helpful, user friendly, and control friendly option. After this gets merged, then either @d-burg or I should make a PR changing the default to D_perp and removing the chi_perp_e and chi_perp_i options. Let me know if you have any questions on this. Thanks!

@logan-nc

logan-nc commented Oct 5, 2026

Copy link
Copy Markdown
Collaborator

At a high level, this looks good. I agree that if someone has chi_phi and chi_perp_e, they likely have chi_perp_i on hand. If they don't have things, or are lazy then letting them assert a "reasonable" tau_E for their device is a great way to get results quick.

@d-burg d-burg left a comment

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks -- the switch is cleanly done, the default path is the same arithmetic as before, and the rename is complete. I checked the new models against the papers: D_perp matches Fitzpatrick (2022) Eqs. 4–6 and 16, c_beta matches Park (2022) and is the exact χ⊥ → 0 limit of the coded D_perp, and P_phi matches Fitzpatrick (2023) Eq. 143. tau_E does not yet match Eq. 131 (inline at line 400). I also compared the old and new example file: only chi_i and the provenance text were added, and every pre-existing dataset is byte-identical.

A mix of top level feedback from myself and claude:

Blocking

  1. Docs CI. Two docstring lines parse as Markdown links (suggestion inline at LayerParameters.jl:220).
  2. Rebase onto develop. This now conflicts in 7 files since #298 landed. Note for the binary file: #298 raised chi_e in TkMkr_D3Dlike_Hmode_kinetic.h5, and this branch's copy predates that, so chi_i needs to be added to develop's file rather than taking this branch's version, otherwise the chi_e change is undone. The chi_i/chi_e ratios in your table will change accordingly.

On making D_perp the default and removing chi_perp_e / chi_perp_i

I agree D_perp is the better-founded model, but I'd settle a few things before switching:

  • D_perp is Fitzpatrick's derived perpendicular particle diffusivity: c_β²η⊥/μ₀ plus (2/3)(1 − c_β²)(n T_e′χ⊥,e + n T_i′χ⊥,i)/p′, built from the electron and ion energy diffusivities. By comparison chi_perp_e is an electron-only proxy, so I'd describe it in the docstring and release note as the backward-compatible choice rather than "the" P⊥. One caveat for the default: where the temperature profiles are flat, Eq. 16 reduces to the c_beta floor, which your table shows moves γ by 60–130 %, so flat-temperature regions of a real profile would change answers a lot.
  • tau_E needs fixing before it is advertised as the quick option: Eq. 131's τ_E is the local r²/χ_E, not the whole-plasma confinement time, so as documented this mode is off by r̂_s²/κ_e (details inline).
  • χ⊥,e stays a required input in every mode because it sets the critical-Δ width, so it can't go away even if the model option does. This PR is already a breaking rename of a key every SLAYER deck sets; I'd rather agree the end state now so users migrate once.
  • Please add one harness case with a non-default P_perp_model (your table already has the reference numbers), so a later default change has a baseline.

Data question. The example file's provenance attribute now names DIII-D 147131 at 2300 ms as the source of chi_i; the previous text named no shot. Has releasing this AOT-derived profile in a public repository been cleared under the DIII-D data policy? Since the provenance is being rewritten anyway, it should also use the dataset names chi_e and chi_phi, not "chi_perp" and "chi_tor".

Small, outside the diff

  • docs/src/kinetic_forces.md documents the kinetic HDF5 schema without chi_i and calls chi_e "the" perpendicular heat diffusivity; the run_slayer docstring likewise mentions only chi_e/chi_phi.
  • An old deck with chi_perp = ... is rejected by the generic unknown-key error, which is the important part, but a one-line hint when chi_perp is among the unknown keys ("renamed to chi_perp_e") would save users a search through several dozen names.
  • Fitzpatrick (2022) and (2023) are cited for three of the six models but aren't in docs/resources; adding them would let later reviewers check the equations.

Comment on lines +220 to +223
- `chi_perp_i` -- ion perpendicular heat diffusivity [m²/s] (used by `P_perp_model=:chi_perp_i`)
- `dlnn`, `dlnTe`, `dlnTi` -- radial log-gradients of n, T_e, T_i in any common radial variable
(only their ratios enter; required by `P_perp_model=:D_perp`)
- `tau_E` -- whole-plasma energy confinement time [s] (required by `P_perp_model=:tau_E`)

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These two lines are what fails the Documentation job: [m²/s] (used by ...) and [s] (required by ...) parse as Markdown links (the parser skips the space between ] and (), and this docstring is spliced into inner_layer.md. A comma after the unit bracket, as the build_slayer_inputs docstring already does, removes both. The suggestion also adds :D_perp to the models that use chi_perp_i.

Suggested change
- `chi_perp_i` -- ion perpendicular heat diffusivity [m²/s] (used by `P_perp_model=:chi_perp_i`)
- `dlnn`, `dlnTe`, `dlnTi` -- radial log-gradients of n, T_e, T_i in any common radial variable
(only their ratios enter; required by `P_perp_model=:D_perp`)
- `tau_E` -- whole-plasma energy confinement time [s] (required by `P_perp_model=:tau_E`)
- `chi_perp_i` -- ion perpendicular heat diffusivity [m²/s], used by `P_perp_model=:chi_perp_i` and `:D_perp`
- `dlnn`, `dlnTe`, `dlnTi` -- radial log-gradients of n, T_e, T_i in any common radial variable
(only their ratios enter; required by `P_perp_model=:D_perp`)
- `tau_E` -- whole-plasma energy confinement time [s], required by `P_perp_model=:tau_E`

Comment on lines +410 to +416
# D_⊥ from Fitzpatrick, Phys. Plasmas 29, 032507 (2022), Eqs. 4-6 and 16 (η_⊥ = η_∥), with
# η_s/(1+η_s) = dlnT_s/(dlnn+dlnT_s) and τ = (T_e/T_i)(1+η_e)/(1+η_i) written flat-density safe.
frac_e = dlnTe / (dlnn + dlnTe)
frac_i = dlnTi / (dlnn + dlnTi)
tau_F = (t_e / t_i) * (dlnn + dlnTe) / (dlnn + dlnTi)
D_perp = c_beta^2 * eta / MU_0 +
(2 / 3) * (1 - c_beta^2) * (frac_e * tau_F / (1 + tau_F) * chi_perp_e + frac_i / (1 + tau_F) * chi_perp_i)

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I checked this against Fitzpatrick (2022) Eqs. 4–6 and 16: the expression and the docstring formula match the paper. Two small things.

The intermediate ratios have a singularity the formula itself does not: with a vanishing ion (or electron) pressure gradient, frac_i and tau_F both go infinite and the product gives Inf·0 = NaN, so this throws although Eq. 16 is finite there. Writing the weights over one common denominator is the same arithmetic (η_e τ_e/(1+η_e) = T_e·dlnTe/dp and likewise for ions), leaves only p′ = 0 singular, and makes the physical content visible: the bracket is (n T_e′ χ⊥,e + n T_i′ χ⊥,i)/p′.

Suggested change
# D_⊥ from Fitzpatrick, Phys. Plasmas 29, 032507 (2022), Eqs. 4-6 and 16 (η_⊥ = η_∥), with
# η_s/(1+η_s) = dlnT_s/(dlnn+dlnT_s) and τ = (T_e/T_i)(1+η_e)/(1+η_i) written flat-density safe.
frac_e = dlnTe / (dlnn + dlnTe)
frac_i = dlnTi / (dlnn + dlnTi)
tau_F = (t_e / t_i) * (dlnn + dlnTe) / (dlnn + dlnTi)
D_perp = c_beta^2 * eta / MU_0 +
(2 / 3) * (1 - c_beta^2) * (frac_e * tau_F / (1 + tau_F) * chi_perp_e + frac_i / (1 + tau_F) * chi_perp_i)
# D_⊥ from Fitzpatrick, Phys. Plasmas 29, 032507 (2022), Eqs. 4-6 and 16 (η_⊥ = η_∥); the weights η_s τ_s/(1+η_s) equal n T_s'/p'.
dp = t_e * (dlnn + dlnTe) + t_i * (dlnn + dlnTi)
D_perp = c_beta^2 * eta / MU_0 +
(2 / 3) * (1 - c_beta^2) * (t_e * dlnTe * chi_perp_e + t_i * dlnTi * chi_perp_i) / dp

It would also help to say in the docstring what Eq. 16 implies at the edges, since users may be surprised: where T_e and T_i are flat, the χ terms vanish and P⊥ falls to the c_beta value (≈ 0.01 here), and where the species pressure gradients have opposite signs (hollow density) the weights can exceed 1. A test at a flat-temperature point would pin the first.

chi_perp::Real, chi_tor::Real,
chi_perp_e::Real, chi_tor::Real,
m::Integer, n::Integer,
chi_perp_i::Real=1.0, tau_E::Real=NaN,

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

tau_E and the log-gradients default to NaN so that a model which needs them fails loudly, but chi_perp_i defaults to 1.0 m²/s, so :chi_perp_i or :D_perp called without it returns a plausible-looking number. I'd default it to NaN here and in build_slayer_inputs, and move the isfinite(...) && ... > 0 check out of the :D_perp branch to just after the if block so it applies to P_perp for every model: today a zero, negative or NaN χ⊥,i (a partly-zero profile, or spline undershoot) gives P⊥ ≤ 0 in :chi_perp_i mode with no error.

+ `:D_perp` -- τ_R/τ_⊥ with τ_⊥ = r_s²/D_⊥ and D_⊥ from Eq. 16 of Fitzpatrick, Phys. Plasmas 29,
032507 (2022), with η_⊥ = η_∥:
D_⊥ = c_β² η/μ₀ + (2/3)(1 − c_β²)[η_e τ_e/(1+η_e) χ⊥,e + η_i τ_i/(1+η_i) χ⊥,i]
+ `:c_beta` -- P_perp = C² with C = c_β (J.-K. Park, Phys. Plasmas 29, 072506 (2022), Eq. 10, with the

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Small citation fix: in Park (2022) the C² term is in Eq. (9) and is defined in the text just below Eq. (10); Eq. (10) itself carries P, which is our P_tor. I'd cite "Eq. 9, with C² = c_β² + (1 − c_β²)K and K = 0" here and at line 423. On line 230, please also define τ_e = τ/(1+τ), τ_i = 1/(1+τ), η_s = dlnT_s/dlnn, and note that this τ is (T_e/T_i)(1+η_e)/(1+η_i), the opposite orientation to the tau = T_i/T_e field of SLAYERParameters.

Comment thread src/Tearing/Runner/run_slayer.jl Outdated
chi_perp_e = loaded.chi_perp_e === nothing ? control.chi_perp_e : loaded.chi_perp_e
chi_tor = loaded.chi_tor === nothing ? control.chi_tor : loaded.chi_tor
(loaded.chi_perp === nothing || loaded.chi_tor === nothing) && @warn(
(loaded.chi_perp_e === nothing || loaded.chi_tor === nothing) && @warn(

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In tau_E mode P_tor is replaced by τ_R/τ_E, so the file's chi_phi and the scalar chi_tor are ignored; that is only stated in the slayer_parameters docstring. Please say it in the SLAYERControl docstring and the TOML annotation too (suggestion on the example deck), and log one @info when the mode is active. This warning should also depend on the model: as written it reports a chi_tor fallback in tau_E mode, and a chi_perp_e fallback in c_beta, P_phi and tau_E mode with dc_type = "none", where those values do not enter the result.

chi_perp_e = 1.0 # fallback only; the kinetic file supplies χ⊥,e(ψ)
chi_perp_i = 1.0 # fallback only; the kinetic file supplies χ⊥,i(ψ)
chi_tor = 1.0 # fallback only; the kinetic file supplies χ_φ(ψ)
P_perp_model = "chi_perp_e" # Perpendicular Prandtl number model: "chi_perp_e", "chi_perp_i", "P_phi", "D_perp", "c_beta", or "tau_E" (needs tau_E [s])

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
P_perp_model = "chi_perp_e" # Perpendicular Prandtl number model: "chi_perp_e", "chi_perp_i", "P_phi", "D_perp", "c_beta", or "tau_E" (needs tau_E [s])
P_perp_model = "chi_perp_e" # Perpendicular Prandtl number model: "chi_perp_e", "chi_perp_i", "P_phi", "D_perp", "c_beta", or "tau_E" (needs tau_E [s]; also replaces the χ_φ-based P_tor)

# η_e = 2, η_i = 3, T_e = T_i ⇒ τ = 3/4, τ_e = 3/7, τ_i = 4/7
d1 = slayer_parameters(; base..., P_perp_model=:D_perp, grads..., chi_perp_i=2.0)
D_hand = d1.c_beta^2 * d1.eta / MU_0 + (2 / 3) * (1 - d1.c_beta^2) * ((2 / 3) * (3 / 7) * 1.0 + (3 / 4) * (4 / 7) * 2.0)
@test d1.P_perp ≈ d1.tau_r * D_hand / d1.rs^2 rtol = 1e-12

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

With T_e = T_i this check passes whichever way round the temperature ratio is written, and the struct's own tau is T_i/T_e, so I'd add a T_e ≠ T_i point (hand numbers: τ = 3/2, τ_e = 3/5, τ_i = 2/5; not run). One build_slayer_inputs test with P_perp_model=:D_perp would also cover the gradient forwarding from the profile splines, which is otherwise untested. I haven't run this test myself so this is just a suggestion, you should confirm on your end whether it checks out.

Suggested change
@test d1.P_perp ≈ d1.tau_r * D_hand / d1.rs^2 rtol = 1e-12
@test d1.P_perp ≈ d1.tau_r * D_hand / d1.rs^2 rtol = 1e-12
# T_e = 2T_i ⇒ τ = 3/2, τ_e = 3/5, τ_i = 2/5: fixes the orientation of the temperature ratio
d2 = slayer_parameters(; base..., t_e=2000.0, P_perp_model=:D_perp, grads..., chi_perp_i=2.0)
D_hand2 = d2.c_beta^2 * d2.eta / MU_0 + (2 / 3) * (1 - d2.c_beta^2) * ((2 / 3) * (3 / 5) * 1.0 + (3 / 4) * (2 / 5) * 2.0)
@test d2.P_perp ≈ d2.tau_r * D_hand2 / d2.rs^2 rtol = 1e-12

Comment thread src/Tearing/Runner/Control.jl Outdated
ctrl.P_perp_model in _VALID_P_PERP_MODELS ||
throw(ArgumentError("SLAYERControl: P_perp_model=$(ctrl.P_perp_model) " *
"not in $(_VALID_P_PERP_MODELS)"))
ctrl.P_perp_model !== :tau_E || (ctrl.tau_E !== nothing && ctrl.tau_E > 0) ||

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Missing, non-positive and NaN tau_E are caught here, but inf gets through and only fails later inside slayer_parameters, after the equilibrium and stability stages. Adding isfinite makes the two checks agree.

Suggested change
ctrl.P_perp_model !== :tau_E || (ctrl.tau_E !== nothing && ctrl.tau_E > 0) ||
ctrl.P_perp_model !== :tau_E || (ctrl.tau_E !== nothing && isfinite(ctrl.tau_E) && ctrl.tau_E > 0) ||

isfinite(tau_E) && tau_E > 0 ||
throw(ArgumentError("slayer_parameters: P_perp_model=:tau_E needs a positive tau_E, got $tau_E"))
# P_⊥ = P_φ = τ_R/τ_E (Fitzpatrick, Phys. Plasmas 30, 092512 (2023), Eq. 131)
P_tor = tau_r / tau_E

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I checked this against Fitzpatrick (2023), Phys. Plasmas 30, 092512. Eq. 131 is P⊥ = P_φ = τ_R/τ_E, but the τ_E there is the local energy confinement time τ_E(r) = r²/χ_E with χ_E uniform (Eqs. 81–83), not the whole-plasma confinement time. Eqs. 140–142 make this explicit: P⊥(r̂) = (τ_R0/τ_E0)·ĵ_oh(r̂) with τ_E0 = a²/χ_E, so the r² cancels and P⊥ = μ₀χ_E/η(r_s). The whole-plasma time from W_e/P_e or a scaling law is τ̄_E = κ_e a²/χ_E (Eqs. 37–38), with κ_e = (3/2)∫₀¹ T̂_e r̂ dr̂ / [−dT̂_e(1)/dr̂].

With tau_E documented as the whole-plasma confinement time, this code's P⊥ differs from the paper's by r̂_s²/κ_e. For a T_e profile like (1 − r̂²), κ_e ≈ 0.19, so the 3/1 and 4/1 surfaces come out about 2–4× larger than Eq. 131 gives, and the error grows outward.

Two ways to make it match: take the whole-plasma τ̄_E and use P⊥ = τ_R(r_s)·κ_e a²/(r_s² τ̄_E), which needs κ_e from the T_e profile and the minor radius; or take a uniform χ_E (= κ_e a²/τ̄_E) as the input, which makes this mode the same as the χ-based P⊥ with one number for every surface. Either way, the docstring should say which confinement time is meant and cite Eqs. 81–83 alongside Eq. 131.

@logan-nc

logan-nc commented Oct 17, 2026 •

Copy link
Copy Markdown
Collaborator

Thanks for the thorough review @d-burg
I think @ebursch will tackle these changes post paternity leave.

I'd just like to reinforce that DIII-D data cannot be included in this open repo directly. The approach for the "DIIID-like" example should be to cite a paper if available and make a roughly matching functional form of the profile needed. Ideally, note the functional form somewhere in annotations / docs so it's clear its not a direct profile fit to data (eg "y = y0 (1 -x^2) where y0=3.0 roughly matches [Name Journal 2026]" or some such).

d-burg and others added 4 commits October 19, 2026 16:42
Brings in the post-solve TearingProblem API, the TOML renames, Z_eff(ψ) from the kinetic file
and the layer-overlap cap. The P_perp_model options now reach build_slayer_inputs through
TearingProblem, and the layer-overlap scan takes the same chi_perp_e, chi_perp_i, tau_E and
P_perp_model as the SLAYER stage, through one helper that also warns only about a missing χ the
selected model uses.

Review changes folded into the resolved hunks:
- D_perp writes the Eq. 16 weights over one denominator, n T_s'/p', so a flat species
  temperature no longer gives 0·∞.
- chi_perp_i has no default (NaN), and every P_perp model now rejects a non-finite or
  non-positive P_perp.
- tau_E is the whole-plasma confinement time: a uniform χ_E = (3/16)a²/τ̄_E sets P_perp and
  P_tor at each surface (Fitzpatrick 2023, Eqs. 37-38, 81-83, 131), with a taken from the
  boundary on the same radial label as r_s. Before, P_perp = τ_R/τ̄_E was too large by
  r̂_s²/κ_e. A non-finite tau_E is caught at validation.
- Docstrings: units no longer parse as Markdown links, Park (2022) is cited at Eq. 9, and
  Eq. 16's symbols and the orientation of its temperature ratio are defined.
- An old deck's chi_perp key gets a hint that it is now chi_perp_e.

The kinetic file is rebuilt from develop's copy, so #298's chi_e raise is kept. chi_e and
chi_phi are replaced, and chi_i added, as analytic profiles in psi_N rather than profiles
derived from a measured discharge; n, T and omega_E are unchanged.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
… and the tau_E conversion, with a D_perp harness case

Adds a T_e ≠ T_i point that fixes the orientation of Eq. 16's temperature ratio, the
flat-temperature c_β² floor, the missing and zero chi_perp_i errors, and the τ̄_E → χ_E
conversion, both in slayer_parameters and through build_slayer_inputs (gradients from the
profile splines, minor radius from the boundary). The new diiid_slayer_n1_d_perp case reruns
the SLAYER deck with P_perp_model = D_perp, so the non-default models have a baseline.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…inetic file and the P_perp models

States the three analytic forms with their sources: Hong et al. (2023) for the magnitude of χ_e
and χ_i, Yoshida et al. (2012) for χ_φ tracking χ_i, and Callen et al. (2010) for the
pedestal barrier the forms do not resolve. Replaces the stability docs' note that separating
the electron and ion channels was deferred with the P_perp_model options, and adds chi_i to the
kinetic schema.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…x label

Every P_perp model is μ₀χ/η, so r_s cancels and the only radial label left in tau_E mode is
the minor radius a in χ_E = (3/16)a²/τ̄_E. The outboard midplane chord, the default rs_method,
is shortened by the Shafranov shift and ignores elongation; the flux label √(2ψ_t/B₀) is the
area-equivalent radius of the boundary and reduces to the minor radius on circular surfaces.
tau_E mode now uses it whatever rs_method is.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
d-burg and others added 3 commits October 19, 2026 17:04
…erp models

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…nd pin the truncation deck to chi_perp_e

[SLAYER] P_perp_model now defaults to D_perp, the perpendicular particle diffusivity of
Fitzpatrick (2022) Eq. 16, which is the coefficient the layer equations and the Eq. 100 width
call for. On the DIII-D-like profiles it is 0.22-0.32 of chi_e at the outer rational surfaces,
so SLAYER growth rates move and the layers are about 30% narrower. A run that selects D_perp or
chi_perp_i without a chi_i profile or a chi_perp_i value now stops with a message naming the
three ways to supply one. slayer_parameters keeps chi_perp_e as its own default, since it has no
profile gradients.

Under D_perp the layers on the DIII-D-like profiles no longer overlap inside psi = 1, so the
truncation deck pins chi_perp_e to keep exercising the layer-overlap cap; it still caps psilim
at 0.99498 with q = 6 the last surface. The D_perp harness case becomes diiid_slayer_n1_chi_perp_e,
a baseline for the electron-channel model.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
@d-burg d-burg changed the title InnerLayer.SLAYER - FEATURE! - Add perpendicular Prandtl number options and an ion heat diffusivity profile InnerLayer.SLAYER - FEATURE! - Build P_perp from the perpendicular particle diffusivity by default and add an ion heat diffusivity profile Oct 20, 2026

This branch has not been deployed

No deployments
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

changed-results Results move or an interface breaks - read before upgrading feature New capability

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants